#include"sm_system.h"
#include"help.h"
#include<fstream>
#include<iostream>
#include<cmath>//abs,sqrt
#include<algorithm>//min,max

using std::ifstream;
using std::abs;
using std::min;
using std::max;

namespace SM
{
  void System_Tz_SM_M::setupOrbitals()
  {
    vector<OrbitalType> OrbitalsTemp;
    vector<int> OrdersTemp;
    //split different m orbitals from System_Tz_SM
    for(int i=0;i<pSystem->Orbitals.size();i++)
      {
	int l=pSystem->Orbitals[i].l;
	int jj=pSystem->Orbitals[i].jj;
	int tz=pSystem->Orbitals[i].tz;
	DataType e=pSystem->Orbitals[i].e;
	for(int mm=-jj;mm<=jj;mm+=2)
	  {
	    OrbitalsTemp.push_back(OrbitalType(l,jj,mm,tz,e));
	    OrdersTemp.push_back(i);
	  }
      }
    //store the Orbitals so that protons and neutrons are separated, proton first.
    //notice that we use this property when generating configurations in Configs.
    //the orbitals with same j but adjacent m is put together, this is important to simplify the cal. of diagonal element of J^2 between configs, see getJ2 in Configs.
    Orbitals.clear();
    Orders.clear();
    Orbitals.reserve( OrbitalsTemp.size() );
    Orders.reserve( OrdersTemp.size() );
    for(int i=0;i<OrbitalsTemp.size();i++)
      {
	if(OrbitalsTemp[i].tz==-1)
	  {
	    Orbitals.push_back(OrbitalsTemp[i]);
	    Orders.push_back(OrdersTemp[i]);
	  }
      }
    maxProtons=Orbitals.size();
    for(int i=0;i<OrbitalsTemp.size();i++)
      {
	if(OrbitalsTemp[i].tz==1)
	  {
	    Orbitals.push_back(OrbitalsTemp[i]);
	    Orders.push_back(OrdersTemp[i]);
	  }
      }
    maxNeutrons=Orbitals.size()-maxProtons;
  }
  
  void System_Tz_SM_M::setupTwoBodyStates()
  {
    TwoBodyStates.clear();
    for(int i=0;i<Orbitals.size();i++)
      for(int j=i+1;j<Orbitals.size();j++)
	{
	  int Tz=(Orbitals[i].tz + Orbitals[j].tz)/2;
	  int Par=(Orbitals[i].par + Orbitals[j].par)%2;
	  int M=(Orbitals[i].mm + Orbitals[j].mm)/2;;
	  TwoBodyStates.push_back(TwoBodyStateType(i,j,Tz,Par,M));
	}
    setupTwoBodyChannels();//noting that the twobody states are sorted after this calling.
  }

  const int  System_Tz_SM_M::MatEledim=System_Tz_SM_M::VSystem::MatEledim;//initialize of static variable
  void System_Tz_SM_M::setupTwoBodyMat()
  {
    setupIndiceIn2BStates();
    int totalChannels=TwoBodyChannels.size();
    TwoBodyMat.resize(totalChannels);
    for(int channel=0;channel<totalChannels;channel++)
      {
	int num2BInChannel=TwoBodyChannels[channel].size();
	TwoBodyMat[channel].clear();
	TwoBodyMat[channel].resize( (num2BInChannel*(num2BInChannel+1))/2 );
	for(int bra=0;bra<num2BInChannel;bra++)
	  for(int ket=bra;ket<num2BInChannel;ket++)
	    {
	      int ab=TwoBodyChannels[channel][bra];
	      int cd=TwoBodyChannels[channel][ket];
	      int a=TwoBodyStates[ab].a;
	      int b=TwoBodyStates[ab].b;
	      int c=TwoBodyStates[cd].a;
	      int d=TwoBodyStates[cd].b;
	      TwoBodyMatEle MatEle(MatEledim);
	      get2BmatOnNeed(a,b,c,d,MatEle);
	      set2Bmat(channel,bra,ket,MatEle);
	    }
      }
  }

  void System_Tz_SM_M::setupIndiceIn2BStates()
  {
    int totalOrbitals=Orbitals.size();
    IndiceIn2BStates.resize((totalOrbitals*(totalOrbitals-1))/2);
    for(int i=0;i<TwoBodyStates.size();i++)
      {
	int a=TwoBodyStates[i].a;
	int b=TwoBodyStates[i].b;
	//a<b
	IndiceIn2BStates[mapab(a,b)]=i;
      }
  }

  void System_Tz_SM_M::get2BmatAt(int a,int b,int c,int d,TwoBodyMatEle & MatEle) const
  {
    MatEle.setZero();
    if(a==b || c==d) return;//pauli principle
    int phase=1;
    //make sure a<b,c<d,(ab)<=(cd)
    if(a>b)
      {
	swap(a,b);
	phase*=-1;
      }
    if(c>d)
      {
	swap(c,d);
	phase*=-1;
      }
    int Mab=(Orbitals[a].mm + Orbitals[b].mm)/2;
    int Mcd=(Orbitals[c].mm + Orbitals[d].mm)/2;

    if( Mab!=Mcd ) return;
    int ab=IndiceIn2BStates[mapab(a,b)];
    int cd=IndiceIn2BStates[mapab(c,d)];
    if(ab==-1||cd==-1) return;
    int channelab=ChannelIndice[ab];
    int channelcd=ChannelIndice[cd];
    //channelab should = channelcd
    if(channelab!=channelcd) return;
    int bra=TwoBodyStateIndice[ab];
    int ket=TwoBodyStateIndice[cd];
    if(bra>ket)
      swap(bra,ket);
    get2Bmat(channelab,bra,ket,MatEle);
    MatEle*=phase;
  }

  void System_Tz_SM_M::get2BmatOnNeed(int a,int b,int c,int d,TwoBodyMatEle & MatEle) const
  {
    MatEle.resize(MatEledim);
    MatEle.setZero();

    if(a==b||c==d) return;//pauli
    int jja=Orbitals[a].jj;
    int mma=Orbitals[a].mm;
    int jjb=Orbitals[b].jj;
    int mmb=Orbitals[b].mm;
    int Tzab=(Orbitals[a].tz+Orbitals[b].tz)/2;;
    int MMab=mma+mmb;
    
    int jjc=Orbitals[c].jj;
    int mmc=Orbitals[c].mm;    
    int jjd=Orbitals[d].jj;
    int mmd=Orbitals[d].mm;    
    int Tzcd=(Orbitals[c].tz+Orbitals[d].tz)/2;
    int MMcd=mmc+mmd;

    int Jmin=max(abs(jja-jjb),abs(jjc-jjd))/2;
    int Jmax=min(jja+jjb,jjc+jjd)/2;

    if(Tzab!=Tzcd) return;
    if(MMab!=MMcd) return;
    if((Orbitals[a].par+Orbitals[b].par+Orbitals[c].par+Orbitals[d].par)%2) return;//parity check
    if(Jmin>Jmax) return;

    int MM=MMab;
    int ordera=Orders[a];
    int orderb=Orders[b];
    int orderc=Orders[c];
    int orderd=Orders[d];
    VSystem::TwoBodyMatEle Temp;
    for(int J=Jmin;J<=Jmax;J++)
      {
	pSystem->get2BmatAt(ordera,orderb,orderc,orderd,J,Temp);
	if(Temp.empty()) continue;
	if(ordera==orderb) Temp*=sqrt(2.);
	if(orderc==orderd) Temp*=sqrt(2.);
	MatEle+=Temp*cg(jja,jjb,J*2,mma,mmb,MM)*cg(jjc,jjd,J*2,mmc,mmd,MM);
      }
  }


  double System_Tz_SM_M::jplus(int a,int b) const
  {
    if(Orders[a]!=Orders[b])
      return 0;
    if( (Orbitals[a].mm-Orbitals[b].mm)/2 ==1 )
      {
	int jj=Orbitals[b].jj;
	int mm=Orbitals[b].mm;
	double j=0.5*jj;
	double m=0.5*mm;
	return sqrt( j*(j+1)-m*(m+1) );
      }
    else
      return 0;
  }

  double System_Tz_SM_M::jminus(int a,int b) const
  {
    if(Orders[a]!=Orders[b])
      return 0;
    if( (Orbitals[b].mm-Orbitals[a].mm)/2 == 1 )
      {
	int jj=Orbitals[b].jj;
	int mm=Orbitals[b].mm;
	double j=0.5*jj;
	double m=0.5*mm;
	return sqrt( j*(j+1)-m*(m-1) );
      }
    else
      return 0;
  }
  
  void System_Tz_SM_M::printOrbitals() const
  {
    for(int i=0;i<Orbitals.size();i++)
      {
	std::cout<<i<<"\t"<<Orbitals[i].l<<"\t"<<Orbitals[i].jj<<"\t"<<Orbitals[i].mm<<"\t"<<Orbitals[i].tz<<"\t"<<Orbitals[i].e<<std::endl;
      }
  }


  void System_Tz_SM_M::printTwoBodyStates() const
  {
    for(int i=0;i<TwoBodyStates.size();i++)
      {
	std::cout<<TwoBodyStates[i].a<<"\t"<<TwoBodyStates[i].b<<"\t"<<TwoBodyStates[i].Tz<<"\t"<<TwoBodyStates[i].Par<<"\t"<<TwoBodyStates[i].M<<std::endl;
      }
  }
  //****************************************************************







  
  void Configs::setup(int MM,int Par,DataType Emax)
  {
    configs.clear();
    //generate proton and neutron separately
    bitset<numBits> start_p,start_n;
    bitset<numBits> last_p,last_n;
    for(int i=0;i<vZ;i++)
      {
	start_p[i]=1;
	last_p[maxProtons-1-i]=1;
      }
    for(int i=0;i<vN;i++)
      {
	start_n[i]=1;
	last_n[maxNeutrons-1-i]=1;
      }
    DataType E0_p,E0_n;
    E0_p=getE(start_p);
    bitset<numBits> config_temp;
    for(int k=0;k<vN;k++)
      config_temp[maxProtons+k]=1;
    E0_n=getE(config_temp);

    int nmax_p=vZ,nmax_n=vN;
    while((Orbitals[nmax_p].e-Orbitals[vZ-1].e)<Emax)
      {
	nmax_p++;
      }
    while((Orbitals[maxProtons+nmax_n].e-Orbitals[maxProtons+vN-1].e)<Emax)
      {
	nmax_n++;
      }
    int maxNum_p=Cnr(nmax_p,vZ);
    int maxNum_n=Cnr(nmax_n,vN);
    
    vector< bitset<numBits> > configs_p;
    vector< bitset<numBits> > configs_n;
    bitset<numBits> config(start_p);
    int count=1;
    while(true)
      {
	if(getE(config) < Emax + E0_p)
	  {
	    configs_p.push_back(config);
	  }
	if(config==last_p || count > maxNum_p)
	  break;
	config=next(config);
	count++;
      }
    config=start_n;
    count=1;
    while(true)
      {
	if(getE(config<<maxProtons)  < Emax + E0_n)
	  {
	    configs_n.push_back(config);
	  }
	if(config==last_n || count> maxNum_n)
	  break;
	config=next(config);
	count++;
      }
    Emax+=E0_p;
    Emax+=E0_n;
    //combine proton and neutron configs togather and choose those the total angular momentum projection is MM/2
    for(int i=0;i<configs_n.size();i++)
      for(int j=0;j<configs_p.size();j++)
	{
	  bitset<numBits> temp(configs_p[j]);
	  for(int k=0;k<maxNeutrons;k++)
	    temp[maxProtons+k]=configs_n[i][k];
	  if(getMM(temp)==MM && getPar(temp)==Par && getE(temp)< Emax)
	    configs.push_back(temp);
	}
    dim=configs.size();
    // for(int i=0;i<configs.size();i++)
    //   std::cout<<i<<"\t"<<configs[i]<<std::endl;
  }

  int Configs::getMM(const bitset<numBits> & config) const
  {
    int MM=0;
    for(int i=0;i<maxProtons+maxNeutrons;i++)
      {
	if(config.test(i))
	  MM+=Orbitals[i].mm;
      }
    return MM;
  }

  int Configs::getPar(const bitset<numBits> & config) const
  {
    int Par=0;
    for(int i=0;i<maxProtons+maxNeutrons;i++)
      {
	if(config.test(i))
	  Par+=Orbitals[i].par;
      }
    return Par%2;
  }

  DataType Configs::getE(const bitset<numBits> & config) const
  {
    DataType E=0;
    for(int i=0;i<maxProtons+maxNeutrons;i++)
      {
	if(config.test(i))
	  E+=Orbitals[i].e;
      }
    return E;
  }




  
  DataType Configs::getME(int i,int j) const
  {
    DataType temp(0);
    bitset<numBits> diff=configs[i]^configs[j];
    bitset<numBits> diffket=diff&configs[j];
    bitset<numBits> diffbra=diff&configs[i];
    int diffCount=diff.count();
    vector<int> pos;
    if(diffCount > 4) return temp;//we only have twobody interaction at most
    else if(diffCount==4)
      {
	findones(diffket,pos,2);
	int c=pos[1];
	int d=pos[0];
	findones(diffbra,pos,2);
	int a=pos[1];
	int b=pos[0];
	temp=get2B(a,b,c,d);
	int sign=count(configs[i],a,b) + count(configs[j],c,d);
	temp=sign%2? -temp:temp;
      }
    else if(diffCount==2)
      {
	findones(diffket,pos,1);
	int b=pos[0];
	findones(diffbra,pos,1);
	int a=pos[0];
	int sign=count(configs[i],a,0) + count(configs[j],b,0);
	temp=get1B(a,b);
	for(int k=0;k<maxProtons+maxNeutrons;k++)
	  {
	    if(configs[i].test(k))
	      temp+=get2B(a,k,b,k);
	  }
	temp=sign%2?-temp:temp;
      }
    else if(diffCount==0)
      {
	for(int k=0;k<maxProtons+maxNeutrons;k++)
	  {
	    if(configs[i].test(k))
	      {
		temp+=get1B(k,k);
		for(int l=0;l<maxProtons+maxNeutrons;l++)
		  {
		    if(configs[i].test(l))
		      {
			temp+=0.5*get2B(k,l,k,l);
		      }
		  }
	      }
	  }
	temp+=get0B();
      }
    return temp;
  }


  double Configs::getJ2(int i,int j) const
  {
    double temp(0);
    bitset<numBits> diff=configs[i]^configs[j];
    bitset<numBits> diffket=diff&configs[j];
    bitset<numBits> diffbra=diff&configs[i];
    int diffCount=diff.count();
    vector<int> pos;
    if(diffCount > 4) return temp;//we only have twobody interaction
    else if(diffCount==4)
      {
	findones(diffket,pos,2);
	int c=pos[1];
	int d=pos[0];
	findones(diffbra,pos,2);
	int a=pos[1];
	int b=pos[0];
	temp=pSystem->jminus(b,d)*pSystem->jplus(a,c)+pSystem->jminus(a,c)*pSystem->jplus(b,d);
	temp-=pSystem->jminus(a,d)*pSystem->jplus(b,c)+pSystem->jminus(b,c)*pSystem->jplus(a,d);
	int sign=count(configs[i],a,b) + count(configs[j],c,d);
	temp=sign%2? -temp:temp;
      }
    else if(diffCount==0)
      {
	//J-J+ matrix element between configs = sum_{hp}(J-)_hp (J+)_ph , p(h) is orbital (not) occupied by this config. because the orbital with same j but with adjacent  m is adjacent, we simplify to use one loop only
	for(int ii=0;ii<maxNeutrons+maxProtons;ii++)
	  {
	    if(configs[i].test(ii))
	      {
		int jj=ii+1;
		if(jj>=maxNeutrons+maxProtons) continue;
		// for(int jj=0;jj<maxNeutrons+maxProtons;jj++)
		//   {
		if(configs[i].test(jj)) continue;
		double j=0.5*Orbitals[ii].jj;
		double m=0.5*Orbitals[ii].mm;
		temp+=j*(j+1)-m*(m+1);
		//		  }
	      }
	  }
	double M=getMM(configs[i])*0.5;
	//Jz^2+Jz
	temp+=M*(M+1);
      }
    return temp;
  }



  
  //****************************************************************


  

}

